function draw_frustum
% define the parameters
P = 1000;
a = 1;
h = 0.1;
E = 1e6;
nu = 0.2;
b = (3*(1-nu*nu)/h/a/h/a)^0.25;
D = E * h * h * h / 12 / (1-nu*nu);

% calculate the exact solution
w = @(x)(-exp(-b.*abs(x))./(2*b*b*b*D).* (P/4) .* (cos(b.*abs(x)) + sin(b.*abs(x))));
x_exact = -3:0.01:3;
w_exact = w(x_exact);

% input the FEM solution
x_FEM = -3:0.1:3;
w_FEM = 10*[0.00000e+00 3.05289e-09 9.85813e-09 2.37457e-08 ...
    4.67763e-08 7.81023e-08 1.11329e-07 1.31027e-07 ...
    1.09189e-07 3.48949e-09 -2.39402e-07 -6.70639e-07 ...
    -1.31025e-06 -2.09704e-06 -2.81608e-06 -2.81608e-06 ...
    -1.91485e-06 1.56829e-06 8.70602e-06 2.05175e-05 ...
    3.69309e-05 5.54230e-05 6.91758e-05 6.50732e-05 ...
    2.23630e-05 -8.64803e-05 -2.90316e-04 -6.06251e-04 ...
    -1.01583e-03 -1.42896e-03 -1.63844e-03 -1.42896e-03 ...
    -1.01583e-03 -6.06251e-04 -2.90316e-04 -8.64803e-05 ...
    2.23630e-05 6.50732e-05 6.91758e-05 5.54230e-05 ...
    3.69309e-05 2.05175e-05 8.70602e-06 1.56829e-06 ...
    -1.91485e-06 -3.01203e-06 -2.81608e-06 -2.09704e-06 ...
    -1.31025e-06 -6.70639e-07 -2.39402e-07 3.48949e-09 ...
    1.09189e-07 1.31027e-07 1.11329e-07 7.81023e-08 ...
    4.67763e-08 2.37457e-08 9.85813e-09 3.05289e-09 ...
    0.00000e+00];

% plot
figure(1)
plot(x_exact,w_exact,'r-',x_FEM,w_FEM,'k--');
legend('exact solution','FEM solution');

figure(2)
subplot(1,2,1)
plot(x_exact(1:100),w_exact(1:100),'r-',x_FEM(1:11),w_FEM(1:11),'k--');
legend('exact solution','FEM solution');
subplot(1,2,2)
plot(x_exact(101:200),w_exact(101:200),'r-',x_FEM(12:21),w_FEM(12:21),'k--');
legend('exact solution','FEM solution');

end